clc,clear;
close all;
warning off;
%数据处理
% 导入数据
filename = 'newdata.xlsx';  
data1 = readtable(filename);
% 提取年份和省份数据
years = data1.x__;
pdata = data1(:, 2:end); 
nyears=length(years);
% 转换表格数据为矩阵
pd = table2array(pdata);
% 输入矩阵 X 
P1 = pd;
P2 = xlsread('data.xlsx','卫生人员');
P2 = P2(:,2:end);
P3 = xlsread('data.xlsx','床位');
P3 = P3(:,2:end);
% 输出矩阵 Y 
Q1=xlsread('data.xlsx','门诊人数');
Q1 = Q1(:,2:end);
Q2=xlsread('data.xlsx','门诊费用');
Q2 = Q2(:,2:end);
Q3=xlsread('data.xlsx','住院人数');
Q3 = Q3(:,2:end);
Q4=xlsread('data.xlsx','住院费用');
Q4 = Q4(:,2:end);
%处理
u=Q1(:,end);
Q1=[u,Q1(:,1:end-1)];
u=Q2(:,end);
Q2=[u,Q2(:,1:end-1)];
u=Q3(:,end);
Q3=[u,Q3(:,1:end-1)];
u=Q4(:,end);
Q4=[u,Q4(:,1:end-1)];
% 环境矩阵 Z
Z1 = xlsread('data.xlsx','人均GPD');
Z1 = Z1(:,2:end);
Z2 = xlsread('data.xlsx','人口总数');
Z2 = Z2(:,2:end);
%对环境矩阵处理
Z1 = (Z1 - mean(Z1)) ./ std(Z1); % 标准化处理
Z2 = (Z2 - mean(Z2)) ./ std(Z2); % 标准化处理

 %% 计算3阶段DEA
for k=1:nyears
        P=[P1(k,:);P2(k,:);P3(k,:)];
        Q=[Q1(k,:);Q2(k,:);Q3(k,:);Q4(k,:)];
        S=[Z1(k,:),Z2(k,:)];
        % 计算 DEA 效率
        %BBC
        n=size(P',1);m=size(P,1);q=size(Q,1);
        w = [];
        for i = 1:n
           f = [zeros(1,n) 1];              % 定义目标函数
           A = [P -P(:,i); -Q zeros(q,1)];  % 指定不等式约束
           b = [zeros(1,m) -Q(:,i)']'; 
           Aeq = [ones(1,n) 0];             % 定义等式约束
           beq = 1;
           LB =[zeros(n+1,1)];              % 指定下界
           UB = [];
           w(:,i) = linprog(f,A,b,Aeq,beq,LB,UB);   % 模型求解 
        end
        thetaa=w(n+1,:)';%评价
        % SFA回归分析
        % 计算输入松弛变量
        opti = thetaa' .* P;  % 最优输入 = 效率值 * 实际输入
        slack = P - opti;  % 输入松弛变量 = 实际输入 - 最优输入
        [m1, n1] = size(slack);  % n1 是决策单元数量，m1 是输入变量数量
        betaes = zeros(size(S, 2) + 2, m1);  % 存储每个输入变量的回归系数
        for j = 1:m1
            % 初始参数估计
            inparam = [zeros(size(S, 2), 1); 1;1];
            
            % 进行SFA回归
            options1 = optimoptions('fminunc', 'Algorithm', 'quasi-newton', 'Display', 'iter', 'SpecifyObjectiveGradient', false);
            [estparam, ~, ~, ~, ~, hessian] = fminunc(@(params) sfa(params, S, slack(:, j)), inparam,options1);
            
            % 存储估计的参数
            betaes(:, j) = estparam;
            %计算标准误差
            cmatrix = inv(hessian);  % 通过逆 Hessian 矩阵计算协方差矩阵
            stder(:,j) = sqrt(diag(cmatrix));  % 标准误差是协方差矩阵的对角线元素的平方根
        end
        % 提取回归系数 beta
        beta = betaes(1:end-2,:);
        % 计算 t 值(每年的要分别查看)
        t = betaes ./stder;
        % 提取随机误差和效率误差的标准差
        sigvest = betaes(end-1, :);
        siguest = betaes(end, :);
        % 调整后的松弛变量
        aslack = zeros(size(slack));
        for j = 1:m
            adjust = S * beta(:, j);
            aslack(:, j) = slack(:, j) - adjust;
        end
        % 调整输入数据 X(P0)
         P0 = P - aslack;
         % 确保调整后的输入数据合理，未出现负值或异常值
         P0 = max(P0, 0);  % 将负值修正为 0
        % 使用调整后的输入数据进行DEA分析
        %CCR
        n=size(P0',1);m=size(P0,1);s=size(Q,1);
        A=[-P0' Q'];
        b=zeros(n,1);
        LB=zeros(m+s,1);UB=[];
        for i=1:n
              f0=[zeros(1,m)  -Q(:,i)'];
              Aeq=[P0(:,i)',zeros(1,s)];beq=1;
              w0(:, i)=linprog(f0,A,b,Aeq,beq,LB,UB);
              E0(i,i)=Q(:,i)'*w0(m+1:m+s,i);
        end
        theta0(k,:)=diag(E0)';%评价(综合)
        %BBC
        n=size(P0',1);m=size(P0,1);q=size(Q,1);
        w = [];
        for i = 1:n
           f = [zeros(1,n) 1];              % 定义目标函数
           A = [P0 -P0(:,i); -Q zeros(q,1)];  % 指定不等式约束
           b = [zeros(1,m) -Q(:,i)']';
           Aeq = [ones(1,n) 0];             % 定义等式约束
           beq = 1;
           LB =[zeros(n+1,1)];              % 指定下界
           UB = [];
           w(:,i) = linprog(f,A,b,Aeq,beq,LB,UB);   % 模型求解 
        end
        theta1(k,:)=w(n+1,:)';%评价(纯)
   
end
theta2=theta0./theta1;
%误差处理
n=size(theta2,1);
m=size(theta2,2);
for i=1:n
    for j=1:m
        if theta2(i,j)>1
            theta2(i,j)=1;
            theta1(i,j)=theta0(i,j);
        end
    end
end
%% 计算 Moran's I 指数
W=readmatrix('经纬度.csv');
W=W(1:end-3,2:3);
W0=W;
W=dist(W);
thetaL=[theta0(:,2:end)]';%除去全国的综合技术
n = size(thetaL, 1); % DMU 数量
Wsum = sum(W(:)); % 权重矩阵 W 的元素之和
Moran_I = zeros(1, size(thetaL, 2)); % 存储每个时间段的 Moran's I 指数
for t = 1:size(thetaL, 2)%至少2个对象
    E_t = thetaL(:, t); % 当前时间段的 DEA 效率值
    E_mean = mean(E_t); % 当前时间段的平均效率值
    
    num = 0; % 公式中的分子
    denom = sum((E_t - E_mean).^2); % 公式中的分母
    for i = 1:n
        for j = 1:n
            num = num + W(i, j) * (E_t(i) - E_mean) * (E_t(j) - E_mean);
        end
    end
   Moran_I(t) = (-n / Wsum) * (num / denom);
end
% 计算每年的标准差
stdev = std(thetaL, 0, 1); 

% 设置年份序列
years = 1990:2023; 

% 绘制图表
figure;
% 绘制全局莫兰指数曲线
yyaxis left;
plot(years, Moran_I, '-o', 'LineWidth', 1.5, 'DisplayName', '全局莫兰指数');
ylabel('全局莫兰指数');
% 绘制标准差曲线
yyaxis right;
plot(years, stdev, '-s', 'LineWidth', 1.5, 'DisplayName', 'DEA效率值标准差');
ylabel('标准差');
% 添加图例
legend('show', 'Location', 'best');
% title('全局莫兰指数与DEA效率值标准差随时间变化');
xlabel('年份');
grid on;
hold off;


%% 局部莫兰指数
%标准化
W= W ./ sum(W, 2);
Moran_I_local = zeros(n, size(thetaL, 2)); % 存储每个时间段的局部莫兰指数
for t = 1:size(thetaL, 2)
    E_t = thetaL(:, t); % 当前时间段的 DEA 效率值
    E_mean = mean(E_t); % 当前时间段的平均效率值

    for i = 1:n
        % 计算局部莫兰指数的分子
        num = 0;
        for j = 1:n
            num = num + W(i, j) * (E_t(j) - E_mean);
        end
        
        % 计算局部莫兰指数
        denom = (E_t(i) - E_mean)^2;
        Moran_I_local(i, t) = (E_t(i) - E_mean) * num / denom;
    end
end

names=readcell('经纬度.csv');
names=names(2:end-3,1);
% 设定要绘制的年份
year_index = 34; % 对应的年份索引

% 获取指定年份的局部莫兰指数
Moran_I_year = Moran_I_local(:, year_index);

% 计算全局均值
gmean = mean(Moran_I_year);

% 划分四个象限
% 象限 1: 高-高 (HH) -> 绿色圆圈
% 象限 2: 高-低 (HL) -> 红色方形
% 象限 3: 低-高 (LH) -> 蓝色三角形
% 象限 4: 低-低 (LL) -> 黑色菱形

% 初始化分类
HH = []; HL = []; LH = []; LL = [];
HH_names = {}; HL_names = {}; LH_names = {}; LL_names = {};

% 遍历每个省份，进行分类
for i = 1:length(Moran_I_year)
    if Moran_I_year(i) >= gmean
        if Moran_I_year(i) >= gmean
            HH = [HH; W0(i,1), W0(i,2)];
            HH_names{end+1} = names{i};
        else
            HL = [HL;W0(i,1), W0(i,2)];
            HL_names{end+1} = names{i};
        end
    else
        if Moran_I_year(i) >= gmean
            LH = [LH;W0(i,1), W0(i,2)];
            LH_names{end+1} = names{i};
        else
            LL = [LL; W0(i,1), W0(i,2)];
            LL_names{end+1} = names{i};
        end
    end
end
% 绘制散点图
figure;
if ~isempty(HH)
    scatter(HH(:,1), HH(:,2), 100, 'g', 'o', 'filled'); 
    hold on;
    text(HH(:,1), HH(:,2), HH_names, 'VerticalAlignment','bottom', 'HorizontalAlignment','right');
end
if ~isempty(HL)
    scatter(HL(:,1), HL(:,2), 100, 'r', 's', 'filled');
    text(HL(:,1), HL(:,2), HL_names, 'VerticalAlignment','bottom', 'HorizontalAlignment','right');
end
if ~isempty(LH)
    scatter(LH(:,1), LH(:,2), 100, 'b', '^', 'filled');
    text(LH(:,1), LH(:,2), LH_names, 'VerticalAlignment','bottom', 'HorizontalAlignment','right');
end
if ~isempty(LL)
    scatter(LL(:,1), LL(:,2), 100, 'k', 'd', 'filled');
    text(LL(:,1), LL(:,2), LL_names, 'VerticalAlignment','bottom', 'HorizontalAlignment','right');
end
% 添加图例
legend({'高-高 (HH)', '高-低 (HL)', '低-高 (LH)', '低-低 (LL)'}, 'Location', 'best');
% 设置地图的标题和标签
% title(['各省份局部莫兰指数分布图 - ', num2str(1990 + year_index - 1), '年']);
xlabel('纬度');
ylabel('经度');
% 显示网格
grid on;
hold off;


%% 基尼系数
regions = [1,1,1,2,2,1,2,2,1,1,1,2,1,...
           2,1,2,2,2,1,1,1,3,...
           3,3,3,3,3,3,3,3,3]; % 1为东部、2为中部、3为西部
theta=theta0(:,2:end);
years = size(theta, 1);  % 获取年份数（34年）
provinces = size(theta, 2);  % 获取省份数（31省）
% 初始化结果存储
G_w = zeros(years, 1); % 总的区域内差异
G_nb = zeros(years, 1); % 总的区域间差异
G_t = zeros(years, 1); % 总的超变密度
G = zeros(years, 1); % 总体基尼系数

% 逐年计算
for year = 1:years
    thetayear = theta(year, :); % 获取某年的theta数据
    % 按照区域划分 theta 到不同的区域矩阵中
    thetaeast = thetayear(regions == 1);
    thetacentral = thetayear(regions == 2);
    thetawest = thetayear(regions == 3);
    theta_combined = [thetaeast'; thetacentral'; thetawest'];
    % 计算全国平均效率
    mtheta = mean(theta_combined);
    % 初始化各个基尼系数部分
    G_wtemp = 0; % 区域内差异
    G_ntemp = 0; % 区域间差异
    G_temp = 0; % 超变密度
    % 计算每个区域的 G_w
    for j = 1:3 % 遍历每个区域
        ij = (regions == j); % 当前区域的索引
        tj = theta_combined(ij); % 当前区域的效率值
        nj = sum(ij); % 当前区域的省份数量
        G_wj = 0;
        for i = 1:nj
            for r = 1:nj
                G_wj = G_wj + abs(tj(i) - tj(r));
            end
        end
        % 将每个区域的 G_wj 归一化，并按照区域权重汇总到总的 G_w
        G_wtemp = G_wtemp + (G_wj / (2 * nj^2)) * (nj / numel(theta_combined));
    end
    % 计算每个区域的 G_nb
    for j = 1:3 % 外层区域
        ij = (regions == j);
        tj = theta_combined(ij);
        mean_j = mean(tj); % 当前区域的均值
        for h = 1:3 % 内层区域
            if j ~= h
                ih = (regions == h);
                th = theta_combined(ih);
                meanh = mean(th); % 另一区域的均值
                nj = sum(ij); % j区域的省份数量
                nh = sum(ih); % h区域的省份数量
                % 区域间差异
                G_ntemp = G_ntemp + abs(mean_j - meanh) * (nj * nh) / (2 * numel(theta_combined)^2);
            end
        end
    end
    % 计算每个区域的 G_t
    for j = 1:3 % 遍历每个区域
        ij = (regions == j);
        tj = theta_combined(ij);
        nj = sum(ij);
        G_tj = 0;
        for i = 1:nj
            G_tj = G_tj + abs(tj(i) - mtheta);
        end
        % 将每个区域的 G_t_j 归一化，并按照区域权重汇总到总的 G_t_total
        G_temp = G_temp + (G_tj / (2 * nj)) * (nj / numel(theta_combined));
    end
    % 计算总体基尼系数 G
    G_w(year) = G_wtemp;
    G_nb(year) = G_ntemp;
    G_t(year) = G_temp;
    G(year) = (G_wtemp + G_ntemp + G_temp);
end
% 计算每个部分的贡献率
cG_w = G_w ./ G * 100; % 区域内差异的贡献率
cG_nb = G_nb ./ G * 100; % 区域间差异的贡献率
cG_t = G_t ./ G * 100; % 超变密度的贡献率

% 绘制图表
figure;
yyaxis left; % 左侧Y轴
% 绘制总体基尼系数曲线
plot((1:years) + 1989, G, '-d', 'DisplayName', '总体基尼系数 G', 'LineWidth', 1.5);
ylabel('总体基尼系数 G');
yyaxis right; % 右侧Y轴
% 绘制贡献率曲线
plot((1:years) + 1989, cG_w, '-o', 'DisplayName', 'G_w 区域内差异', 'LineWidth', 1.5);
hold on;
plot((1:years) + 1989, cG_nb, '-s', 'DisplayName', 'G_nb 区域间差异', 'LineWidth', 1.5);
plot((1:years) + 1989, cG_t, '-^', 'DisplayName', 'G_t 超变密度', 'LineWidth', 1.5);
ylabel('贡献率 (%)');
% 添加图例
legend('show', 'Location', 'best');
% 添加标题和标签
xlabel('年份');
% title('总体基尼系数及其组成部分的贡献率随时间变化');
grid on;
hold off;

%% 区域间差异(基尼系数)
% 初始化存储每个区域内差异的数组
G_weast = zeros(years, 1);   % 东部区域内差异
G_wcl = zeros(years, 1);% 中部区域内差异
G_wwest = zeros(years, 1);   % 西部区域内差异

% 逐年计算
for year = 1:years
    thetayear = theta(year, :); % 获取某年的theta数据
    
    % 按照区域划分 theta 到不同的区域矩阵中
    thetaeast = thetayear(regions == 1);
    thetacentral = thetayear(regions == 2);
    thetawest = thetayear(regions == 3);
    
    % 计算东部区域内差异
    G_wjeast = 0;
    njeast = numel(thetaeast);
    for i = 1:njeast
        for r = 1:njeast
            G_wjeast = G_wjeast + abs(thetaeast(i) - thetaeast(r));
        end
    end
    G_weast(year) = G_wjeast / (2 * njeast^2);
    
    % 计算中部区域内差异
    G_wjcentral = 0;
    njcentral = numel(thetacentral);
    for i = 1:njcentral
        for r = 1:njcentral
            G_wjcentral = G_wjcentral + abs(thetacentral(i) - thetacentral(r));
        end
    end
    G_wcl(year) = G_wjcentral / (2 * njcentral^2);
    
    % 计算西部区域内差异
    G_wjwest = 0;
    nj_west = numel(thetawest);
    for i = 1:nj_west
        for r = 1:nj_west
            G_wjwest = G_wjwest + abs(thetawest(i) - thetawest(r));
        end
    end
    G_wwest(year) = G_wjwest / (2 * nj_west^2);
end

% 绘制区域内差异图
figure;
plot((1:years) + 1989, G_weast, '-o', 'DisplayName', '东部区域内差异 G_weast', 'LineWidth', 1.5);
hold on;
plot((1:years) + 1989, G_wcl, '-s', 'DisplayName', '中部区域内差异 G_wcentral', 'LineWidth', 1.5);
plot((1:years) + 1989, G_wwest, '-^', 'DisplayName', '西部区域内差异 G_wwest', 'LineWidth', 1.5);

% 添加图例
legend('show', 'Location', 'best');

% 添加标题和标签
xlabel('年份');
ylabel('区域内差异');
% title('东部、中部、西部区域内差异随时间变化');
grid on;
hold off;

%% Kernel 密度估计
eMatrix =theta0(:,2:end);
% 创建年份和效率值的坐标
years = repmat((1:34)', 1, 31);
eValues = eMatrix(:);
years = years(:);
% 生成网格数据
[effGrid, yearGrid] = meshgrid(linspace(0, 1, 100), 1:34);
% 对每个网格点进行插值估算密度
densityGrid = zeros(size(effGrid));
for i = 1:34
    % 提取当前年份的数据
    eValues_year = eMatrix(i, :);
    % 核密度估计
    [f, xi] = ksdensity(eValues_year, 'Kernel', 'epanechnikov', 'Bandwidth', obandwidth(eValues_year));
    % 移除重复的xi值
    [xi_unique, unique_idx] = unique(xi);
    f_unique = f(unique_idx);
    % 在网格上插值
    densityGrid(i, :) = interp1(xi_unique, f_unique, effGrid(i, :), 'linear', 0);
end
% 绘制三维曲面图，X轴为效率值，Y轴为年份
figure;
surf(effGrid, yearGrid + 1989, densityGrid, 'EdgeColor', 'none');  % 曲面图，无边缘线
colormap(jet);  % 使用jet颜色映射
colorbar;
% title('医疗资源配置效率动态演进');
xlabel('医疗资源配置效率');
ylabel('年份');
zlabel('核密度');
set(gca, 'XLim', [0 1]);  % 限制X轴（效率值）在0到1的范围
set(gca, 'YTick', 1990:5:2025);  % 设置Y轴为离散的年份刻度
grid on;
% 设置视角以更好地展示三维结构
view([-30, 20]);

%% 调用函数
%求省份之间坐标
function distance=dist(x)
    n=length(x);
    distance=zeros(n,n);
    for i=1:n
        for j=1:n
            if i==j
                distance(i,j)=0;
            else
                lat1=x(i,1)*pi/180;
                lat2=x(j,1)*pi/180;
                lon1=x(i,2)*pi/180;
                lon2=x(j,2)*pi/180;
                dlon=lon2-lon1;
                dlat=lat2-lat1;
                a=sin(dlat/2)^2+cos(lat1)*cos(lat2)*sin(dlon/2)^2;
                c=2*atan2(sqrt(a),sqrt(1-a));
                R=6371; %地球半径
                distance(i,j)=R*c;
            end
        end
    end
end

% 带宽选择的辅助函数
function bw = obandwidth(data)
    n = numel(data);
    sigma = std(data);
    bw = 1.06 * sigma * n^(-1/5);
end

function negLogL = sfa(params, Z, slack)
    beta = params(1:end-2);  % 环境变量回归系数
    sigma_v = params(end-1); % 随机误差项的标准差
    sigma_u = params(end);   % 效率误差项的标准差

    % 计算残差
    residual = slack - Z * beta;

    % 计算对数似然函数
    sigma = sqrt(sigma_v^2 + sigma_u^2);
    lambda = sigma_u / sigma_v;
    eps1=1;
    % 计算对数似然
    logPhi = log(normcdf(lambda * residual / sigma)+eps1);
    logf = -0.5 * log(2 * pi) - log(sigma) - 0.5 * (residual / sigma).^2 + logPhi;
    negLogL = -sum(logf);  % 负对数似然
end


% 
